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Turbulent mixing of multiphase flow 

By Y.-N. Young 

J. Ferziger, F. E. Ham and M. Herrmann 


1. Motivation and objectives 

Mixing of multi-phase fluids is common in diverse research fields, and understanding 
of such an important phenomenon is useful to a range of engineering applications such as 
the polymer blender or control and designs of fluid processing. In this paper we adopt the 
phase-field modeling approach to investigate capillary induced effects on mixing. Phase- 
field modeling has been applied to simulation of multi-phase flow (Chella & Vinals 1996; 
Jasnow & Vinals 1996; Jacqmin 1999) due to its attractive aspect of easy numerical 
implementation. More recently it has also been applied to simulations of solidification of 
dendritic alloys (Zhao et al. 2003) and polymer blenders (Roths et al. 2002). 

On macroscopic scales, the interface between two fluids is infinitely sharp. In phase- 
field modeling, the fluid properties are modeled (assumed) to change smoothly over a 
small layer across the interface on mesoscopic scales (Chella & Vinals 1996; Jasnow & 
Vinals 1996). Such a diffuse interface of separation is first considered in a liquid near or at 
its critical point, where the fluid property within the boundary needs to be modeled and 
resolved (Siggia 1979). However, phase-field modeling is not restricted to liquids at or near 
the critical point (Chella & Vinals 1996; Jasnow & Vinals 1996; Jacqmin 1999; Young 
et al. 2003). In fact, one can certainly envision that on the time scale of macroscopic 
fluid motions, the thermodynamic potential may depend on the gradient of the order 
parameter, which designates the fluid phase (or density). Accordingly, the macroscopic 
jump in the normal stress across the interface due to surface tension can be modeled 
as a mesoscopic continuum body force proportional to gradient of surface potential. 
On the mescoscopic scales, such a continuum surface tension force varies smoothly in 
the transition region around the interface. Naturally such a smooth continuum model 
of fluid interface brings great numerical convenience. Recently, phase-field modeling is 
made exactly energy-conserving by adding a thermodynamical constraint (Jamet et al. 
2002). We did not include this constraint in this work as the correction is small for our 
cases of interest. 

An interesting finding in multiphase fluids stirred by chaotic mixing flow is the self- 
similar drop size distribution (Berthier et al. 2001; Muzzio et al. 1991 b). After rescaling 
the drop size distribution with respect to the characteristic size, all the drop size distri- 
butions collapse to a universal functional form for various combinations of parameters, 
such as viscosity, capillary number, and stirring flow strength. Some explanation to the 
self-similarity is provided in (Muzzio et al. 1991a, c), and it is not obvious if similar 
collapse of size distribution can also be found in a more general turbulent flow, where 
the distribution of stretching and stretching rate is randomly distributed. In the case of 
stochastic mixing (Lacasta et al. 1995), the existence of a characteristic domain size is 
an indication that similar collapse of drop size distribution may be found. However, no 
report on the self-similar drop size distribution is given in Lacasta et al. (1995). Thus 
we conduct numerical simulations of multiphase fluids stirred by two-dimensional turbu- 
lence to assess the possibility of self-similar drop size distribution in turbulence. In our 
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turbulence simulations, we also explore the non-diffusive limit, where molecular mobility 
for the interface is vanishing. Special care is needed to transport the non-diffusive inter- 
face. Numerically, we use the particle level set method to evolve the interface. Instead of 
using the usual methods to calculate the surface tension force from the level set function, 
we reconstruct the interface based on phase-field modeling, and calculate the continuum 
surface tension forcing from the reconstructed interface. 

This paper is organized as follows. In §2 we formulate the coupled Navier-Stokes-Cahn- 
Hilliarcl system and present results from direct numerical simulations of the coupled 
system stirred by turbulent flows. In §3.1 we formulate the multiphase fluid problem 
combining the particle level set method (for interface tracking) with phase-field modeling 
(for surface tension force). Numerics for the particle level set method are summarized in 
§3.2, and validation is provided in §3.3. In §3.4 we present results of drop dynamics from 
direct numerical simulations, and finally we provide some conclusion in §4. 


2. Turbulent mixing in the Navier-Stokes-Cahn-Hilliard system 

We generalize the system in Berthier et al. (2001) to include the effect of surface 
tension on stirring flows. The phase- field function C denotes the phase of the fluid: 
C = ±1 corresponds to fluid phase 1 and 2, respectively. C is described by the Cahn- 
Hilliard equation, which is coupled to the incompressible, two-dimensional Navier-Stokes 
equations in our formulation. 

— + u • VC = /3MV 2 V>, (2.1) 

— + u • Vu = -VP + zA7 2 u - CVV> + F + F d , (2.2) 

V • u = 0, (2.3) 

where 

^ C) ~ ~ 1 ? 2C ' ir ( C ') = i( C ' 2 - 1 ) 2 - (2- 4 ) 

In equation 2.2, F is a random forcing at wavenumber k in the range n\ < fc/fco < ri 2 , 
where fco = 27 t/1 and l is the dimension of the computation domain. = — Au (A = 0.1 
in all our simulations) is the drag force dissipating energy at large scales. 

The positive, constant coefficient A4 is the molecular mobility, v is the kinematic 
viscosity, a and (3 are related to the interfacial properties: The interfacial thickness £ ~ 
a/P, and the surface tension coefficient a = ^po^/a/3 = \J a 0ao, where po is the 
uniform fluid density and op is the surface tension force when \J aft = 1. We can non- 
dimensionalize the above equations with respect to a characteristic length scale l (of the 
order of interfacial thickness), a characteristic velocity t’o , and the corresponding time 
1/vq- The dimensionless equations are 


— + u • VC = OV 2 [C 3 - C - V 2 C] , (2.5) 

yj L 

r\ -i 

+ u • Vu = -VP + — V 2 u - [C 3 - C - V 2 C] + F' + F' d , (2.6) 

V • u = 0. (2.7) 


The dimensionless coefficient D = (3A4 /vqI is the inverse Peclet number, Re = Vq l/v is 
the Reynolds number, and 7 = \faf3jv qu is the capillary number. If the scaling length 
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Figure 1. Energy spectrum of turbulent velocity (solid line), and the k~ 5 ^ 3 inverse cascade 
scaling is the dashed line, fco = = 1/4. 


scale l is the interfacial thickness (/ = f) and the characteristic velocity is the “diffusion 
velocity” across the interface (no = (3M/f), then the Reynolds number Re = f3M/v and 
the capillary number 7 = a/l3A4v. The time scale for diffusion across the interface is 
tc = i 2 /PM. The characteristic hydrodynamic time scale is rs = 1/|S|, where |S| is 
the velocity strain rate. It is convenient to control the interfacial properties using both 
a and (3, thus in the numerical code the dimensional equations 2. 1-2.4 are implemented 
instead. 

For the following simulation results, the domain size of the double-periodic computa- 
tion box is 87 r x 87 t, and the numerical resolution is 512 2 . The small-scale forcing range 
is fixed at wavenumbers 5 < k/ko < 15. At the beginning of each simulation, an initially 
circular interface of radius 1.57T is placed at the center of a two-dimensional turbulent 
flow, and the volume average of the concentration is (C) = 0.41 > 0. The energy spec- 
trum of the turbulent flow is the inverse cascade with a slope of —5/3 as shown in figure 1. 
Throughout the simulations the energy spectrum is little affected by the surface tension 
force from the interface. 

A numerical challenge in turbulent simulations of NS-CH system is that the diffusion 
coefficient has to be small to avoid phase homogenization. If the diffusion coefficient is 
too large, or if the surface tension is too small, drops disappear as the order parameter 
C homogenizes and settles to the phase that is closest to the mean value. Thus, for 
a given turbulent strength and numerical resolution, we optimize the number of drops 
by minimizing the diffusion with enough surface tension force so we can collect drops of 
various sizes. We integrate the turbulent flow over a long duration so the number of drops 
fluctuates many times around the mean value. The area distributions collected from three 
simulations with resolution 512 2 , after rescaled with respect to the average drop area, are 
shown in figure 2(a). The three rescaled area distributions collapse reasonably well in the 
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Figure 2. Size distribution for simulations of turbulent NS-CH flows. Panel (a): Diamonds are 
for a = 2, (3 = 200 and M = 10 -6 . Crosses are for a = 2, (3 = 200 and M = 5 x 10~ ' . Asterisks 
are for a = 3, (3 = 300 and A4 = 10 -6 . The numerical resolution is 512 2 and v = 0.1 for all three 
simulations, and the curve is ~ exp— S/S*. Panel (b): Diamonds are for a = 2, (3 = 200 and 
M = 2 x 10 ' . Asterisks are for a = 1.5, (3 = 150 and M = 2 x 1CP 7 . The numerical resolution 
is 1024 2 and v = 0.1 for all both cases, and the curve is ~ exp —S/S*. 


range 0.3 < S/S* < 20, despite the poor statistics at both ends of the distribution. The 
diffusion coefficients (see caption in figure 2) that we used in the numerical simulations 
are close to the minimum for this numerical resolution (512 2 ). 1024 2 grid points are used 
for two simulations as we decrease the mobility to A4 = 2 x ICC ' , and similar scaling in 
the distributions is also found as shown in figure 2(b). Despite the limited reliable range 
in both figure 2, the general trend suggests that similar scaling may also exist in the 
tubulent mixing case. 


3. Phase-field modeling with zero mobility M. = 0 

In reality the mobility A4 for multiphase flow can be very small, and the equation 
for the phase field concentration C reduces to the advection transport equation. This 
distinguished limit corresponds to the non-diffusive, immiscible two-fluid flow, where the 
interfacial thickness is ideally zero and thus requires special numerical treatment in the 
simulations. 

In the case of zero diffusivity, we use the phase-field modeling as a means to reconstruct 
the interface knowing its location. We use particle level set method to evolve the interface 
between the two phases. In the level set framework, the zero of the scalar level set function 
cf> implies the location of the interface. To retain the structure of the interface, we assume 
that the Cahn-Hilliard potential ip governs the interfacial energetics, and reconstruct the 
interface for a given interfacial thickness £(= \J a./ (3 as in §2). In §3.1 we formulate this 
system in terms of level set representation of the interface with phase-field modeling for 
interfacial structure with zero molecular mobility. In §3.2 we briefly discuss the numerics 
and some validation is given in §3.3. We then present results from direct numerical 
simulations of multiphase mixing in two-dimensional turbulence in §3.4. 

3.1. Formulation 

Basically the equations are exactly the same as those in §2 except the Cahn-Hilliard 
equation is now reduced to an advection equation for C : 





Mixing of multiphase flow 


243 


— + (u-V)C' = 0. (3.1) 

In the case of vanishing molecular mobility across the interface, the advection transport 
equation 3.1 for concentration C is a numerical challenge, and more complication arises 
in the evaluation of the surface tension force. Our remedy is that, despite the non- 
diffusiveness of the phase field, we assume that the interfacial thickness is finite based on 
the free energy if. Naturally, if we adopt the Cahn-Hilliard surface free energy for if, we 
expect the concentration C to behave as C ~ tanh(^//3/a|r — ?’o|), where |r — r 0 | is the 
distance between the point r and the interface (located at ro). A great advantage of this 
approach is that the surface tension force can be calculated directly from the phase-field 
modeling as a continuum force — CV0, which peaks at the interface. 

Thus instead of solving equation 3.1 over the whole computation domain for C, we 
solve the same equation for a level set function 0 only at the interface (where 0 = 0), and 
re-initialize the level set function 0 to a signed distance function away from 0 = 0. We 
replace equation 3.1 with the level set equation valid only at the fluid interface 0 = 0: 

0*0 + u • V0 = 0|^_ 0 . (3.2) 

Away from 0 = 0, inside a band of width ~ 6 dx around 0 = 0, we re-initialize 0 to a 
signed distance function. To achieve this we have to re-initialize the level set function at 
each time step by solving the following equation to a stationary state 


0 T 0 + s<?n(0o)(|V0| - 1) = 0, (3.3) 

where sgn((p 0 ) is the sign of 0(x, t = t 0 ) defined as 


sgn(0o) 


00 

o + |V0 o | 2 da: 2 ’ 


(3.4) 


for a given grid spacing dx. The normal vector and the mean curvature can be calculated 
in terms of the level set function 0: n = V0/|V0| and n = V • n. The concentration 
phase-field C is then reconstructed from the level set function (£, as previously defined, 
is the interfacial thickness): 


C = tanh(|) if |0| < a£, (3.5) 

C = j^j- if |0| > af, 

where a is a constant in the range 2 < a < 5 so that the concentration C transitions 
smoothly from |C| < 1 to |Cj = 1 at the edge of the band around the interface. From 
our fully-resolved simulations of drop dynamics presented here, we conclude that the 
volume integral of C remains constant at all time provided that the level set is accurately 
capturing the interface. 

Numerically there are various ways to calculate the surface tension force from the level 
set function. Physically the surface tension force leads to a pressure jump across the 
interface. In our one-fluid formulation the surface tension is formulated as a body force 
F ct (= — CV0 in the phase field modeling). Ideally, in the limit of zero interfacial thick- 
ness, the surface tension force is non-zero only at the interface: F CT = (ann + d s as)S(<j)), 
where d s a is the Marangoni force, with s the arclength and s the unit vectors along the 
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tangential direction. On a discrete numerical grid, however, the surface tension force F CT 
has to be smoothed over several grid spacings to avoid un-controllable numerical oscil- 
lation. The usual cosine-delta function used to calculate F CT on a uniform grid (Sussman 
et al. 1994) leads to artificial broadening of the interface, which in turn leads to large 
unphysical parasitic currents around the numerically smoothed interface. In our phase 
field modeling, we assume an interfacial structure based on the the Cahn-Hilliarcl surface 
free energy. A continuum surface tension force is then derived based on the surface free 
energy. As shown in (Jacqmin 1999), this continuum surface tension force gives very small 
parasitic currents compared with that from the cosine-delta construction of the interface. 
Recently a numerical scheme has been developed to minimize the parasitic currents (En- 
right et al. 2002). In our approach, the parasitic currents are small but not as minimized 
as in Enright et al. (2002). The turbulent flow in our resolved simulations of drop dy- 
namics is much larger than the parasitic currents, and thus we expect the constraint in 
Enright et al. (2002) to have little effects on the dynamics and drop area distributions 
presented in §3.4. The formulation we proposed for reconstructing the interface from the 
phase-field modeling through level set function (equations 3. 2-3. 5) is somehow similar 
to the splitting of particle advection from particle diffusion in the random walk particle 
method (Ghoniem & Sherman 1985). In random walk particle methods, particles can be 
first transported by the velocity field at their positions, and then perform random walks 
to account for diffusion. Such a numerical scheme is commonly used in particle tracking 
when the particle diffusion is very small. Similarly in our formulation we first evolve 
the interface by solving the advection equation only near the interface using the hybrid 
particle level set method. Once the location of the interface is advanced and the level set 
is re-initialized as a signed distance function, we construct the phase-field C based on 
the steady state solution to equation 2.1 without advecting flow, 

C(C 2 - 1) - -X/ 2 C = 0. (3.6) 

P 

Assuming that the interfacial structure is only a function of the signed distance to the 
interface, we then express C in terms of the level set function <f> in equation 3.5. 

In summary, the system of iso-density, iso-viscosity fluid with an interface separating 
the two phases is described by equations 2.2, 2.3, and 3.2 with the surface tension force 
—CVip, where V’ is the Cahn-Hilliard free energy and the concentration C is calculated 
from equation 3.5. 


3.2. Numerics 

The level set equation (3.2) is solved using the hybrid particle level set method (Enright 
et al. 2002). We use WENO 5th order scheme to calculate the flux in equation 3.2, 
and the TVD-3rd order RK scheme to advance the level set in time. Following Enright 
et al. (2002), Lagrangian particles are added inside a narrow band (T p ) around the 
interface and are transported using the same RK 3rd order scheme. Before and after 
re-initialization (equation 3.3), particle correction (for detail see Enright et al. (2002)) is 
conducted to ensure that the zero level is not shifted numerically. This hybrid level set 
method is at most volume-conserving as long as particle density along the zero level is 
sufficient. In all our simulations, we use at least 4 2 particles for each cell within the (r p ) 
band. Re-seeding of particles is conducted if the average particle density is found to be 
below 75% of the initial particle density. We re-initialize the level set within a second 
band (T r ) enclosing the particle band (r p ). We refer readers to Enright et al. (2002) 
for detail of the algorithm. Standard tests have been performed (Young et al. 2002) and 
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results are identical to those in the original paper (Enright et al. 2002). To couple the 
level set equation(s) to the fluid solver, we advance both the level set <f> and the fluid 
velocity u using the 3-rcl order explicit scheme. The diffusion part in the momentum 
equation is incorporated as integrating factors in spectral space. The surface tension 
force is calculated from the interface at the previous time step, and we use the same 
velocity fields from the previous time step to advect the velocity, interface, and particles. 

3.3. Linear capillary waves 

To demonstrate that the surface tension force is properly captured in the phase-field mod- 
eling, we compute the oscillatory frequency of the capillary waves for a two-dimensional 
drop. The initial condition for </> is the signed distance function perturbed by a sinusoid 
with a small amplitude and a wavenumber n 

t = 0) = r — ro(l + O.O2sin(?r0)), (3.7) 

where r 0 is the drop radius, r = |x — x 0 | is the distance between x and the drop center x 0 , 
and 9 = tan~ 1 ((y — yo)/(x — Xo))- Analytically (following Rayleigh (1892)), the oscillation 
frequency of capillary waves on the interface separating fluids of the same density and 
zero viscosity (1 /Re = 0) can be expressed in terms of surface tension coefficient <r, drop 
radius ?’o, and wavenumber n of the perturbation 



We obtain the linear solutions from simulating the full system (described in section 3.1) 
with small non-linearity: starting from the initial condition u = 0 and ft given in equation 
3.7, the non-linearity remains small at all time because the small disturbance is damped 
by viscosity. In all our simulations in this section, the viscosity v is fixed at v = 0.1, and 
the perturbation amplitude decays exponentially at rate A with a capillary oscillation 
of frequency u>. We compute the capillary frequency to numerically by extracting the 
oscillation from the decaying amplitude. 

We first calibrate the relationship between <7 and (a, (3). From §2 the surface tension 
a is proportional to \faf3. With the ratio (3/ a fixed (the interface thickness is thus fixed 
at ~ 8 grid spacings), we expect the surface tension a to be linearly proportional to 
a (or P). Table 3.3 shows the computed oscillation frequency from simulations and the 
corresponding surface tension coefficient a from equation 3.8. We find that the surface 
tension a (computed from lo in our case via equation 3.8) indeed varies linearly with 

For 1 < n < 6 we compute the oscillation frequency ui from the simulations with 
a = 0.9817 and /3 = 89.1875 (a = 94.32cro from table 3.3). As shown in figure 3, excellent 
agreement is found between the computed and analytical values of u> for all 1 < n < 6. In 
addition we also remark that the parasitic currents are as small as 10 -6 with a numerical 
resolution 128 2 in a square box of size 27 t x 27 t. 

3.4. Drop dynamics in two-dimensional turbulence 
Drop dynamics and the drop break-up in simple flow configurations has been extensively 
investigated both numerically (Zaleski et al. 1995; Cristini et al. 1998; Tauber et al. 2002) 
and experimentally (Muzzio et al. 1991 &) . Results from experiments (Muzzio et al. 1991 6) 
and modeling (Berthier et al. 2001) on passive chaotic mixing of immiscible fluids show 
that the statistics of drop size collapse due to the self-similar breakup and coalescence 
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Perturbation Wavenumber n 


Figure 3. Capillary oscillation frequency w as a function of n. Symbols are simulation data, 
and line is from analytical result for a = 94.32cto and r o = n. 


Table 1 . Parameters a and (3 used in the linear simulations, and the capillary frequency u> com- 
puted from the time-varying perturbation amplitude, n = 2, v = 0.1, and numerical resolution 
is 128 2 . 



a 

P 

U 

cr/<7o 

Case 1 

0.98 

89.18 

3.02 

94.32 

Case 2 

1.96 

178.37 

4.24 

188.64 

Case 3 

2.94 

267.56 

5.22 

282.96 

Case 4 

3.92 

356.75 

6.04 

377.28 


processes. However, there are only a few recent direct numerical simulations on drop 
dynamics in turbulent flows (Tryggvason et al. 2001). Due to the difficulty in tracking the 
interface accurately while conserving fluid mass, it is a great numerical task to capture the 
statistics of drop size distribution after a series of breakups and coalescence in numerical 
simulations. 

Due to numerical diffusion, it is in general a great numerical challenge to conserve 
total volume for each fluid using level set methods. The hybrid particle level set method 
conserves volume almost perfectly in standard tests (Enright et al. 2002) , and is the only 
level set method that can both track the interface accurately and minimize the mass 
loss in a controllable fashion. The particle level set method, in combination with the 
phase field modeling of the surface tension force, is a good tool for tracking interface 
accurately with minimum mass loss. Provided that enough Lagrangian particles are used 
in tracking the interface, the mass loss can be reduced to a tolerable level. Thus we are 
able to accurately evolve the dynamics of drops in a turbulent flow, and our main goal is 
to examine the statistics collected from DNS dataset. In our simulations, less than 1% of 
maximum volume loss is found for all cases presented here. For the following simulation 
results, the combination of a and /? is chosen so that the interface thickness is fixed 
(~ 6 dx at least). The double-periodic computation domain is of size 87t x 87t, and the 
numerical resolution is 512 2 . 
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t=0.6 t=1.9 t=13.1 



Figure 4. Drop deformation and break-up in two-dimensional turbulence, a = 2.95, 
/3 = 118.92, v = 0.2 and A = 0.1. (a) t = 0.6 (b) t = 1.9 and (c) t = 13.1. 


We fixed the range of small-scale forcing at wavenumbers 5 < fc/fc 0 < 15, and the 
large scale friction = — O.lu. At the beginning of each simulation, an initially circular 
interface of radius tt is placed at the center of a two-dimensional turbulent flow. The 
energy spectrum of the turbulent flow is the inverse cascade with a slope of —5/3 as 
shown in figure 1. Throughout the simulations the energy spectrum is little affected by 
the surface tension force from the interface. The well-known effect of bubbles/drops on 
the turbulent energy spectra (Mazzitelli et al. 2003) is also absent here because there is 
no buoyancy/density contrast in our simulations. 

In figure 4 we demonstrate how the initial drop gets deformed, stretched, and finally 
breaks up into smaller drops: In figure 4(a) (t = 0.6) the interface is highly distorted, and 
in figure 4(b) (t = 1.9) several break-ups lead to smaller children drops. At late times 
(t > 10) the drops are small and the surface tension force is sufficient to overcome the 
stretching to prevent further break-ups, and a statistical equilibrium is reached as shown 
in figure 4(c) ( t = 13.1), where there are a lot of small drops of irregular shapes and 
various sizes. 

We can easily calculate the total length of the interface from the level set as 


where 


Lengtli{0 = 0} 


\X7 H (cf>(x.))\dxdy 


6(<M X ))I y<f\dxdy, 




1 

2 

d 

dip 


1 + 


H (</>), 


1 
7 r 


sin( 



(3.9) 


(3.10) 

(3.11) 


with e = 3 ~ 4 dx as in Sussman et al. (1994). Panel (a) in figure 5 shows the evolution 
of the arclength (scaled to the initial length), and panel (b) is the time- history of kinetic 
energy. We find that the early growth in the length is accompanied by a similar growth of 
the total number of drops as several break-up processes have occurred during (0 < t < 3). 
After t ~ 10 (about 2.5 eddy turnover times) the amplification in length saturates and 
fluctuates throughout the simulations. We integrate to at least 15 eddy turnover times 
to attain statistically equilibrium states. The interfacial thickness is fixed at ~ QAdx for 
all three curves in figures 5. With the same small scale random forcing F, large scale 
dissipation F^ and viscous dissipation v = 0.2, the only difference between the three 
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Figure 5. (a) Evolution of total arclength scaled to the initial value, (b) Corresponding 

evolution of number of drops. 


Table 2. Averaged number of drops and arclengths for the three simulations. 



a 

P 

cr/<To 

L/Lo 

N 

Case 1 

1.47 

59.45 

9.35 

7.51 

79.98 

Case 2 

2.94 

119.91 

18.71 

5.50 

39.69 

Case 3 

5.89 

239.83 

37.42 

3.66 

15.18 


simulations is the surface tension coefficient: er = 9.36<to, a = 18.71<to, and a = 37.43cto 
for curves 1, 2, and 3, respectively. (As defined earlier erg is the surface tension coefficient 
when yfaj3 = 1.) The balance between stabilizing surface tension force and stirring 
turbulent flow results in a critical length scale (Lacasta et al. 1995) 


A 


2 

c 


o 

7 r a 


(3.12) 


where u is the turbulence characteristic velocity. Since the total area (volume) of each 
fluid phase has to conserve in two (three)-dimensions, we can estimate the dependence of 
average number of drops ( N ) on the surface tension. For a given characteristic turbulence 
velocity, in two-dimensions, 


1VA 2 ~ total area of fluid phase 1 = constant — >7V~cr _1 . (3.13) 

Similarly in three-dimensions, 

iVA 3 ~ total volume of fluid phase 1 = constant — > N ~ er~ 3//2 . (3.14) 


From equations 3.13 and 3.14, the total circumference (area or arclength in two and three- 
dimensions, respectively) is inversely proportional to the square root of surface tension: 
N A c in two-dimensions and N A 2 in three-dimensions both lead to ~ 1 /y/a. These scaling 
results in two-dimensions are consistent with the time-averaged values of our simulation 
data in figures 5. In table 3.4 we list the time-averaged N and total arclength for all 
three simulations. 

In the statistically equilibrium state the number of drops fluctuates around the mean 
value, we collect enough statistics for the drop area distribution by integrating over a 
long period of time. Following Berthier et al. (2001) and Muzzio et al. (19916), we rescale 
the drop area distribution by the mean area (calculate by taking moments of the area 
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Figure 6. Size distribution for the three cases. Diamonds are for case 1, asterisks are for case 2 
and crosses are for case 3. Open circles are for a different forcing function but same parameters 
as in case 1. The dashed line is ~ exp— (S/S*) . 


distribution functions as in Muzzio et al. (1991 &)). Since the effective diffusion coefficient 
is assumed to be zero in our turbulence simulations, we cannot calculate the mean area 
based on the scaling of the characteristic domain size L* and D* as in Berthier et al. 
(2001). The rescaled area distributions collapse as shown in figure 6. One additional 
simulation data set (open circles) with a different interfacial thickness and turbulent 
random forcing is added in figure 6. For this additional simulation, the surface tension 
a = 18.71(7o and the interfacial thickness £ = 5.0dx and the turbulent random forcing is 
at wavenumbers 5 < k/ko < 7.5, narrower than the other three simulations. However, it 
seems that the collapse works reasonably well regardless of the differences between the 
simulations. 


4. Conclusion 

We have extended the phase-field modeling to the non-diffusive case. Surface free 
energy from the phase-field modeling is used to calculate the surface tension force based 
on the reconstructed interface using the level set function. Numerically the non-diffusive, 
infinitely sharp interface has to be smoothed over a few grid points. Some interfacial 
structure is always assumed in the numerics, and the resultant surface tension force 
inherits the modeled structures. One popular way to reconstruct interface in level set 
method is the cosine-delta function. However, no physical meaning can be related to 
such a construction, and as a result, the parasitic currents are undesirably large. We 
replace the cosine-delta function with the smooth interfacial structure from the stationary 
solution in phase-field modeling, and we calculate the surface tension force from the 
corresponding surface free energy. The only assumption we make in this approach is 
that the reconstructed interface is a function of the signed distance to the interface. We 
assume that “some” molecular dynamics acts on such a fast time scale that the interfacial 
thickness remains the same, even when the flow is vigorously acting on and around 
the interface. The collapse of drop area distribution in our two-dimensional turbulence 
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simulation confirms our conjecture that the self-similar drop area distribution exists for 
more general velocity fields. 

The particle level set method, combined with phase-field modeling, has great promise 
for large scale, three-dimensional numerical simulations of multiphase flow. One of our 
future goals is to generalize our simulations of two-dimensional drop dynamics to three 
dimensions. We also plan to include density and viscosity contrast for multi- fluid using 
the phase-field modeling. Finally, we are also investigating how the phase-field modeling 
can be applied in simulating micro-fluidic multiphase flows carried out in experiments 
(Truesdell et al. 2003). 
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